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ABSTRACT 


This thesis examines the differences of deterministic and stochastic 


LANCHESTER-type combat models. Using an example of square-law attrition, 
solution methods and solutions are described. A new analytic solution 
for equal attrition rate coefficients is given. The numerical comparison 
includes hypotheses about the expected force levels and the variability 


in the expected force levels as a function of time, initial force levels, 


and breakpoint force levels. 
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I. INTRODUCTION 


Combat models are widely used as decision aids in the defense- 
planning process, at least within the NATO alliance. Current operation- 
al combat models are very complex because combat is a very complex 
process. Unfortunately it is difficult (if not impossible) for the 
beginner to understand the modelling approaches, concepts and motivation, 
that may have been used to build such operational models. However, one 
frequently considers a simple model as a paradigm for the development and 
understanding of such complex models. This basic approach will be used 
in this thesis to explore certain issues in the on-going debate about 
the relative merits of stochastic and deterministic combat models. 

A simple model is examined to explore differences between a determin- 
istic and a stochastic approach to a certain type of analytical combat 
model. As already mentioned, combat is a very complex process, but it 
is also a complex random process, which can be supported by many examples 
from military history. Analytical models are abstractions and very often 
simplifications of reality. It seems to be a legitimate question to ask, 
what effects the further abstraction of neglecting the randomness in 
combat may have. At this moment, it should be pointed out that within 
existing operational analytical models, both stochastic and deterministic 
models are used. 

Previous work done by SPRINGALL [9] and CLARK [4] evolved around 
theoretical aspects. Their main concern was to give exact analytical 


solutions and their proofs. CRAIG [5] started to explore the differences 





between stochastic and deterministic models more from the numerical 
point of view, which will be continued in this thesis. 

In the next chapter, a deterministic and stochastic version of a 
differential combat model will be described. The deterministic versions 
are well-known as LANCHESTER's equations of modern warfare, which were 
developed in 1914. Combat models, which model attrition from enemy 
action through a system of differential equations, are usually referred 


to as LANCHESTER-type models of warfare. 





TI. THE PARADIGM MODELS 


A. THE DETERMINISTIC MODEL 

First, LANCHESTER's equations of modern warfare will be briefly 
reviewed and some simple extensions given. 

In 1914 LANCHESTER [7] hypothesized that under "modern conditions" 
in a combat between two homogeneous forces the firepower of the surviving 
weapons of one side can be concentrated on the surviving targets of the 
other side, so that each side's casualty rate is proportional to the num- 


ber of enemy firers. This can be described by the following equations: 


dx _ 

dt — -ay (2.1) 
ay 

dt bx (220) 


with initial conditions 


X(0) = x (2.3) 


0 
Y(O) = Yo (2.4) 


where a is the attrition rate with which the Y-force attrits the 
X-force, similarly for b. X, and Y) are the initial force levels and 
X(t) and Y(t) are the force levels at time t. The force levels, as a 


function of time t, can be written as 


X(t)=x,cosh(Vab t)-Va/b yysinh(Yab t) (2.5) 
¥(t)=y,cosh (ab t)-1b/a x9sinh( Yab t) (2.6) 


The state equation relating initial force levels with force levels at 


some time t can be derived by dividing (2.1) by (2.2), separating 





variables and integrating to yield 


b(xp=X(t)*)=alyg--¥(t)*). (2.7) 


This form of the state equation explains why this model is also referred 
to as the "“square-law" attrition process. WEISS [11] has given a set of 
assumptions under which LANCHESTER's equations for modern warfare may 
apply: 

Al.) Two homogeneous forces are engaged in combat. Every unit 

On a particular side has the same capabilities. The attrition 


rate may be different for the two forces. 


A2.) Each unit on one side is within weapon range of all units 


on the other side. 


A3.) The effects of successive rounds on the target are 


independent. 


A4.) Each unit has perfect knowledge of target locations and 
fires only at live target (one at a time) killing them at a 
constant rate, which does not depend on the number of targets 


alive. 
A5.) Fire is uniformly distributed over surviving targets. 


The above model implies a fight until one force is annihilated. 
Therefore the model will be slightly changed by introducing the concept 


of unit breakpoints, Xbp and V5 » Which are simply force levels at which 


p 
the side, who reaches it first "breaks off" the engagement, leaving the 
other side as a winner. Also, to be more precise, it should be noted 


that negative force levels for breakpoints equal zero or force levels 





less than nonzero breakpoints are impossible. So the deterministic 


LANCHESTER-type combat model with "square-law" attrition takes the 


following form: 


p 
Ypp< V(t) (2.8) 


0 otherwise 


-bx Xpp<X(t) 


Ypps V(t) (8) 


0 otherwise 


with initial conditions 


X(0) = XQ 


¥(0) = Yo 


(25110) 
(2.11) 


The model in this form, equations (2.8) through (2.11), was used for 


comparisons throughout the thesis. 


B. THE STOCHASTIC MODEL 


There are 
type models. 

* The 

* The 

weakening the 

* The 

* The 


Over time may 


several ways to include random variations in LANCHESTER- 
These are: 

attrition rate coefficients may be random variables. 
enemy's initial force level may be a random variable, 
assumption of perfect knowledge. 

breakpoints may be random variables. 

casualty rate is fixed, but the occurrence of casualties 


be random. 


The only random variation considered here will be the random occur- 


rence of casualties over time. Another specification was to choose a 


10 





model similar to the "square-law"attrition in order to allow comparisons. 
In other words, the question to be asked is "how do random fluctuations 
in the occurrence of casualties modify the deterministic results of the 
square-law attrition process?" 

The approach used here was a continuous parameter MARKOV chain model, 
where the time t varies .continuously and the number of combatants on 
each side is a non-negative integer. Let M(t) be the size of the X-force 
at time t with a particular state value m. Let N(t) be the size of the 
Y-force at time t with a particular state value n. Let Mo and Ng be the 


initial force levels and Mp? Np be the breakpoint force levels of X 


p 
and Y respectively. Fig. 1 shows the state space of this MARKOV chain 
model. Note that at a given time t, any state is described by the two 
force levels of the X and Y force. As each side loses units due to 


attrition and no replacements are allowed, it is easy to understand why 


BILLARD [1] referred to this type of process as a "bivariate death process." 


1] 
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For the description of the stochastic square-law attrition process 
corresponding to the two deterministic differential equations (2.8) and 
(2.9), a system of many differential equations, depending on the battle 
termination model, is required. This system will be given for a fixed- 
force-level-breakpoint battle with square-law attrition. 

The following assumptions yield the stochastic square-law attrition 
process. 

Al.) The attrition process depends only on the current system 

State and time and not on the past history (this assumption is 


usually referred to as markovian property). 


a5 one x casualty during the | E 
A2.) The probability no. POAC ieeo bth = ah 


7 one y casualty during the a 
A3.) The probability P Ga ie ee bh 


A4.) The probability of more than one casualty occurring in 
the time interval t to tth is of the order of magnitude o(h), 


where Tim, 90h) /h = 0. 


A5.) No more casualties can occur once m = Mp orn = Np: 


p 


Making the time interval h infinitesimally small, the following 
set of forward CHAPMAN-KOLMOGOROV equations can be developed. Let 
P({t,m,n) be the probability that the system is in state (m,n) at a 
time t. For convenience each equation is related to a region in the 


State space shown in Fig. 2. 


Ts 





For M=Mo and N=No>» Region I 


9 


d = 
qe toMy Mg) = -(anptbmy) P(t,my.no) (25114) 


for m <M<M and Nn=No> Region II 


| 


om 


QO. 
“oO 


(t,m,no) = any P(t,m+1,n,)-(an +bm) P(t m,n 


0 (23) 


0) 


QO. 


ie 


mor n <n<ny and M=Mq » Region III 


Oo 


QO. 
“9 


(t,m,.n) = bm, P(t,my .n+1)-(antbm,) P(t myn) (2.14) 


QO. 


iE 


for m,_<m<m, and n <n<No > Region IV 


D Q bp 


Oo 


QO. 
wo 


(t.sm,n) = anP(t,m+1,n)+bmP(t,m,nt+1)-(antbm)P(t,m,n) 
OZ a) 


O. 


iE 


for m=m), and Nyy <N<No >» Region V 


QO. 
=e) 


(t,m))>n) = anP(t,m, +1 5n) (2.16) 


OQ. 


0 


for n=n and M5 <m<mg » Region VI 


5) a es ee eas 


O. 
9 


qe tsmon ) = bmP(t,m,n, +1) (2 ile) 


bp 


for m=m), and N=Nh 5» Region VII 
fee OPC 


P(t ,m =p otr. dll t (2.18) 


Ape lbp? 


because of the definition of a breakpoint force level. The 
initial condition is 


P(O,m = ]..0 (2.19) 


92p) 


14 








Figure 2 - REGIONS IN THE STATE SPACE 





As P(t,m,n) is a joint probability distribution, the follow- 
ing must also be true 


femP{(t,m,n)21.0 (C220) 


and 


ma P(t,m,n) = 1.0 (ap) 
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IIT. SOLVING THE DIFFERENTIAL EQUATIONS 


A. SOLUTIONS FOR THE DETERMINISTIC MODEL 

Force levels as a function of time were already given in equations 
(2.5) and (2.6). It is relatively easy to obtain analytical solutions 
for such simple deterministic models as described before. On the other 
hand, it is necessary to point out that for models with any degree of 
operational realism, analytical methods for solving the differential 
equations are usually not available. Therefore, some numerical method 
with a digital computer is usually used. TAYLOR [10] has summarized in 
his Appendix C the most widely used numerical methods, a discussion of 


which seems unnecessary at this point. 


B. SOLUTIONS TO THE STOCHASTIC MODEL 

Even for this relatively simple stochastic model with fixed force 
level breakpoints, which are usually nonzero, a complete set of general 
analytical solutions for the differential equations (2.12) through (2.19) 
has not been found. On the other side, given some minor restrictions 
like a fight to the finish or equal attrition rate coefficients, solu- 
tions, or at least solution methods have been proposed which will be 
briefly discussed in the next section. 

First, the method for getting the state probabilities used here wil] 
be described. Numerical solutions were obtained using the fourth-order 
RUNGE-KUTTA method, which is probably one of the best known finite dif- 


ference approximations to ordinary differential equations (next to the 


Wy 





EULER-CAUCHY-method). To increase the accuracy of the overall solutions 
analytical results for certain regions of the state space were substituted. 
These analytical solutions will be stated now. For region I, i.e. no 


casualties on either side, the solution to (2.12) is 
P(t,my M9) = exp -(anytbm,) t (3.1) 


which can be derived by the standard method used for this kind of dif- 
ferential equation. For the boundary cases, region II and III, i.e. 
one of the two sides has not yet had a casualty, TAYLOR [10] has given 
the analytical expression as 


for ME p< msm and n=Ng; Region II 


Meee OO Pt 
P(t,msny) = J faro/oceP*-1}) exp(-(bm tang) t) , (3.2) 
where J = My-m 


fOr Ny <n<ng and M=Mq » Region III 


P(t,mosn) = kr fonoratet*-0 fs exp(-(bmptang) t) , (3.3) 


where K = Ng-n. 

These two equations (3.2) and (3.3) were obtained by 
recursively solving equations (2.12), (2.13) and (2.14) "from 
the top down." 

C. OTHER ANALYTICAL SOLUTIONS 

The solutions or solution methods for getting the state 

probabilities will only be stated for the square-law attri- 


tion process. Only two were used for the numerical work for 


mms thesis. 
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Apparently, one of the "oldest" analytic solutions was 
given by BROWN [3] in 1955 for the general stochastic 
LANCHESTER-type combat model with time independent attrition 
rates. His approach and solution will be briefly outlined 
for square-law attrition. Consider a path from state 
( My » No ) to some state (m,n). This path can be described 
as a sequence of J = ({ Mo -m) zeros and K = ( No -n) ones, 
where a zero corresponds to an X casualty and a one to a Y 
casualty. Using the binary representation of a positive 
integer, one can define to each realization of a battle path 
an integer k given by 

k =d d d 


Rome Kd tK 


th casualty along a battle path 


where dy oy = ] if the r 
corresponding to k is a Y casualty and dy a O otherwise. 
Also let I; K be the set of all positive integers whose 


binary representation contains exactly K ones and J zeros. 


Then 
m = m,-rt x d, . 
Stag 0 j=l Ker 
Yr 
Ne or = Ny - 2 dy og: 


Then BROWN [3] has shown that 


fe 





oT al Ke 


kel) -o yr=(Q 





Bie.m.,n) = 


exp(-iut)-exp(-tl(m, .,n, .)) 


: du oe) 
] Tipe ya oll -1uU 


where i = >| Se 


T(m,n) = antbm, 


? 7 Jk rt] 
kyr T-iu/ (my ron =) 


and 


Ik rt] de rt 1a pth lay ng JBM 

There was no indication that nonzero breakpoints were 
excluded. A discussion of this solution follows in the next 
section. 

About 14 years later, in 1969, CLARK [4] proposed another 
approach which TAYLOR [10] called a “hybrid analytical-numer- 
ical" method. The restriction is that the breakpoints have 
to be zero, i.e. it is a fight to the finish. Although 
proposed for a general time independent attrition function, 
this approach will be outlined for "square-law" attrition. 
Then according to CLARK [4] the state probabilities are given 


by 
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for O<msmy and O<nsnp 


M9 n 
P(t,m,n) = 7 z oa -exp (-(ak+bj) t} , Cas 


j=m ken 
Zor O<msm, and n=0 
gue Tie 
P(t,m,0) = Coen + 2 : ra -exp(-(ak+bj)t), (3.6) 
: j=m e 9 
for m=0 and O<nsny 
0; 40 
mie,0,n) = enol “ee x -exp (-(ak+bj)t) .(3.7) 
9 im = 9 


mmgeac jast there is to remember that P(t,0,0) = 0 for all 
times t. The constants Goan are determined by a system of 


partial difference equations. 


For O<m< jmp and O<n<kSnq 


cm,ntl 
= an Teck bm isk 


a(n-k)+b(m-j) STG!) 


2] 





for O<m<J<m, and O<n= Kno 


m,n mt+i yn 
Ge =anl. 
J sn J sn 
b(m-j) e329) 


for O<m=Jj<my and O<n<k<np 


m,n m,nt] 
an k = bm net 
a(n-k) ees 0) 


for O<m=J<mp and O<n=k<np 


but (m,n) # (my,NQ) 


m m m 
0 m,n m,n 


0 
Pep ms: -=— ba (3.11) 


j=m ke=nt] oP k j=m+] : 


m,n 
with ,” a He) iy, 
0° 0 
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Also for O<msjsmp and O=n<ksnp 


cms | 


Se 


ak+bj 


m,0O 
C5 5k 


0 


Similarly for O=m<jsmo and O<nskésn 


Then for Tsmsmp 
m n 
m,0 0 a m,0 
oo 7 . Gk 
j=m k=] 


Weim ee 
J O,n 
i a 
| j=] ken 


Sealez)) 


(Seal) 


(3.14) 


Cees, 


Though having the publishing date of 1979, the next 


approach was published in June 1980 by BILLARD [1]. 


She 


considered the LANCHESTER-type square-law attrition combat 


model as a pure death-process and applied SEVERO's [8] 


recursive theorem for solving differential equations. As 


before, only a fight to the finish has been considered. 


és 





The first step is to identify each point (m,n) in the 
State space by a counting coordinate k, where 


k = ( mot1)( n,t1) 


0 0 
-m( Nytl)-n : oo) 
Then 
P({t,m,n]= P[(t,k] (eo ln) 


and the differential equations (2.12) through (2.18) take on 
a slightly different form. As an example, (2.15) will be 
given by 
CPt yk) = anP[t,k-n,-1]+bmP[t,k-1] 
dt b ] 3 0 9 
-(antbm)P[t,k]. (3.18) 
The whole set of differential equations was then 


expressed in matrix terms as 


qeh(t) = BP(t) (3.19) 


with a solution given as 


P(t) = Ce (t) (3.20) 
where e(t) is the (mo+1) (ngtt) x1 -vector with elements 
exp(b, t) with b, being the , th diagonal element of the 


Matrix B. The matrix B can be partitioned into submatrices, 
whose m-coordinate is common, due to the ordering defined by 
the counting coordinate k (Equation 3.16). 

Then 


B = (b 


Daas U,V = 1,2...mg+1 
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where the submatrices Diy have the elements 
Diy = (Og MB als be DP>qQ = T,2...ngtl. 
So the matrix B has the elements 


by lP>P) = -a(no-ptl)-b(my-utl ) 
for u = T,2...m and p = 1,2...n 


> P- 1) = b(my-utl) Pole Ul SS. ae Mot 
Dutaond. . On, Nyt] 
and 
by u-7 (PP) = a(ng-p+) for u = 2,3...mg+l 
P = 1,2...ngt1. 


All other elements are zero. 

Thus, the matrix B has at most 3 nonzero entries per row 
or column. The matrix C can be partitioned in the same way. 
Then using SEVERO‘s [8] theorem and the special form of the 
matrix B, only a part of the C-matrix needs to be determined. 
This part will be omitted here, but the final result will 


be given by 


k 
P(t,k) = D_ clk, j)-exp(b,t) Fea 


J21 


where c(k,j) is the (k,j)t" 


element of the solution matrix C. 
The previous two approaches have required that the force 
level breakpoints be zero. Now, a result will be given whose 
restriction is that the attrition rate coefficients be equal, 
but nonzero breakpoints are allowed. For further reference 


it will be called the Equal-Attrition-Rate-Coefficient- 


Solution (EARCS). 


Zo 





Let a = b = f. 


For urs <mSmq and ny <nsnp 


Cimon 


P(t,m,n) = 
Mytny-m-n) | 


eee 0 MNT nxn (-ffmen)t) (3.22) 


for m = Mp and Ny p<MSNq 


See EEE SE ee 


_ fn C(m,_.+1,n) 
P(tsmi 9 sn) = bp : 


J! 
1-exp(-ft(m +14n+k)) 
J k fd b 
° > (-1) 6 f My pt ltn+k 9 (25280) 
k=0 
where J = My *Ng-M, y-t-n 
for n= Mbp and My p<msmg 
ee itt 1) 
P(t,msn).) a bp 
K l-exp (-ft(m+1+n +j)) 
gis : (3.24) 
(-1)d FUT, 5) , , 


7=0 
where K = My Ng-M-n, WT. 


The coefficients C({m,n) satisfy for My p<M<mM and NE <N<ng 


p 
the following partial difference equation 


C(m,n) = nC(m+1,n) + mC(m,nt1) (3.25) 
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with the boundary conditions 
m-m 
(ny) Fife My p< MSM 


C(m,nq) 


and 
Nyon 
om. jn) = (my) TOY Nb 


<ns 
0 ait G) 


0 


This result has been developed using a method verbally pro- 
posed by TAYLOR. The method will now be outlined. ISBELL 
and MARLOW [6] described a stochastic LANCHESTER-type attri- 
tion process with a different attrition function. Instead of 
the attrition of one force being proportional to the number of 
of enemies of the opposing force as in the square-law attri- 
tion (e.g. for the attrition of the M-force let 

A(m,n) = an 
be the attrition rate and 

B(m,n) = bm 
the attrition rate for the N-force with square-law attrition), 
their attrition rates looked like 

A(m,n) = antcm 
and 

een) = bmtdn 
with the restriction that 

ease = btd. 
But with c=d=0 and a=b=f we are back to square-law attrition. 
amis) leads to equations (3.22) and (3.25). Equations (3.23) 
and (3.24) were derived in the following way (e.g. for 3.23). 


Solving equation (3.22) for m = Mp7! yields 


Zi 





( a Ft(m, )#1+n)) 
p ees tl 3n = “exp(- m +]+n 
bp Mytng-m,,-T-n bp 
+ = eee 
Wrest cyto “On bp! (G28) 


substituting for the second factor its BINOMIAL expansion 


J 
J 
eet) 2 ST (yk (") afk (4 )d-k 


K=0 
with J = MytNg-M, 9 | and multiplying through by the third 
factor. Then using equation (2.16), the differential equation 
for m = Mp and Ny p<nsng » and substituting equation (3.26) 
into the extended form, it can now easily be integrated to 


yield equation (3.23). 


tee DISCUSSION 

In the discussion of the analytical solutions outlined in 
the last section, there is one important point. BROWN [3] 
himself points out that unless m is close to Mo and n is 
close to No his result (equation (3.4)) is of “little practi- 
cal interest." Most of the analytical solutions, especially for 
more general LANCHESTER-type models, have little more than 
“symbolic” character. BROWN's solution is a good example 
ooetnat. 

In comparing the solutions given by BILLARD [2] and 
CLARK [4], this author has the feeling that both solutions are 


equivalent and only the representation is different. This 


28 





intuitive guess needs verification. It may be concidence 
that CLARK [4] and SEVERO [8] published their work in the 
Same year. 

The last presented solution (equation (3.22) through 
(3.24)) seems to be relatively handy for use on a computer. 
It has a big advantage over numerical solution methods other 
than its accuracy, because it is an exact result. Like 
CLARK's method, the coefficients have to be calculated only 
Once for a given set of input data. Then, to get the state 
probabilities for a certain point in time you have to make 
only one set of calculations, as opposed to the numerical 
methods where one has to go from time t = 0 to time t = t 


in small time steps and then have only an approximate result. 


SeeeemPORTANCE OF THE STATE PROBABILITIES 

The state probabilities as a function of time are the 
key to calculating several quantities of interest. These 
are expected force levels as a function of time, variances 
and standard deviations in the force levels and also the 
probability of winning. These quantities are necessary to 
legitimately compare the stochastic with the deterministic 
results. 

To get at least a feeling of how the state probabilities 
evolve over time, the joint probability distribution will be 
presented in a 3-D-picture. It is indeed surprising that 


more use has not been made of computer graphics to 
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investigate the dynamics of a stochastic LANCHESTER-type 

combat model. Table 1 gives the data used for the next five 
figures. These figures may be thought of as "Snapshots" of 
the joint probability for the survivors in this battle taken 


at a sequence of increasing times. 
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TABLE 1 


Data for the Numerical Example | 


Force Levels My = 40 Mbp = 0 
Ng = 40 Mbp = 0 
Attrition Rates a = 0.008 M casualties per 


minute and N firer 


b = 0.004 N casualties per 
minute and M firer 
At Times ty =s0 025 te 
to = 0.250 te 
t. = 0.500 te 
ty = 0.750 te 
te = 1.000 te 


where te = 155.81 minutes is the time a 
deterministic battle with the same force 
levels and unit breakpoints ends, i.e. 


x(t) = 0.0 and y(te) = 28.28 
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Figure 3 - FLOT OF JOINT PROBABILITIES P(t,m,n) 


with data according to Table 1 
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rigure & - rLOT OF JOINT PROBABILITIES P(t,m,n,) 
with data according to ‘rable 1 
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rigure 6 = rLOT OF JOINT PROBABILITIES P(t,m,n) 


with data according to Table | 
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P(t,m,n) 
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rigure 7 - PLOT OF JOINT PROBABILITLES P(t,m,n) 


with data according to ‘lable 1 
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At the time t = 0 all probability is located at 
(My »Ng) in the state space (Region I). As time passes, the 
probability mass is distributed over more states, with the 
mode moving away from the diagonal towards the winning side. 
All points in the state space with breakpoints, i.e. (M5520) 
and (mn) for all m and n, are absorbing states, probability 
mass 1s absorbed in that states. The sum of probability mass 
in Region V (see Fig. 2) represents the probability, that the 
N-force wins at that given time, in Region VI that the M-force 
wins. 

The next five figures show a similar sequence of plots for 
the joint probability P(t,m,n) for the force levels M(t) and 
N(t). The data is explained in Table 2. Note the small dif- 
ferences because of the nonzero force level breakpoints. 
Probability mass having reached the breakpoint "piles" up 
there. The state space is reduced by the fixed force level 
breakpoints, but the probability distribution evolves basical- 


ly in the same qualitative manner as in the previous example. 
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TABLE 2 


Data for the Numerical Example 2 


Force Levels Mo = 40 Mp = 8 
No = 40 Nb = 8 
Attrition Rates a = 0.008 M casualties per 
minute and N firer 
b = 0.004 N casualties per 
minute and M firer 
At Times ty =" Oli02 5 te 
t, = 0.250 te 
t = 0.500 te 
ty = 0.750 te 
te = 1.000 te 


where te = 120.68 minutes is the time a 
deterministic battle with the same force 
levels and unit breakpoints ends, i.e. 


x(t) = 8.0 and y(te) = 28.83 
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Figure 9 = PLOT OF JOINT PROBABILITIES P(t,m,n) 


with data according to rable 2 
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rigure 1u - PLOT OF JUINT PROBABILITIES P(t,m,n) 


with data according to ttable 2 
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Figure 11 = PLOT OF JOINT PROBABILITI£S P(t,m,n) 
with data according to Table 2 
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rigure 12 - PLOT OF JOINT PROBAMILITIES P(t,m,n) 


with data according to ‘rable 2 
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IV. EVOLUTION OF THE FORCE LEVELS 


A. THE DETERMINISTIC MODEL 

The evolution of force levels as a function of time, X(t) and Y(t), 
has already been stated in equations (2.5) and (2.6). The next two 
figures show the force levels with different breakpoints. It is easy to 
realize that introducing a nonzero breakpoint does not change the under- 
lying function. In other words, using X(t) as an example, for both 
figures the same curve was used but at the point where the X-force reaches 
its breakpoint, the curve is "cut." From that point in time, there are 
no more changes in the force levels. So introducting a nonzero break- 
point only shifts the discontinuity (marked by DX1 in Fig. 13) up along 
the curve to DX2. 

The probability for one side to win is either one or zero, because 
it is a deterministic model. To easily determine which side is going to 
win, a victory prediction condition can be obtained by solving each force 
level equation (2.5) and (2.6) for the time to reach its breakpoint tXhp 
by substituting X(t)=xp, and tYb5 by substituting V(t)=yp,- Then X will 
win if tYpp < tXbp» which leads to the prediction condition. X will win 
a fixed force breakpoint battle if and only if 
a(yé : Yop) 


b( x4 ~ po) ; (45) 


Xo/Yo a 


This shows that given the initial data one can predict the outcome 
of the battle in terms of force levels and time until the battle 


finishes. 


4a 
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Y(t) 


DX1 


time t 


Figure 13 - FORCE LEVELS OVER TIME 


Deterministic Model with Breakpoint Xp o> O 
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Figure 14 - rORCE LEVELS OVER TIME 


Deterministic Model with Breakpoint Xn? 0 
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B. THE STOCHASTIC MODEL 
Every possible state (m,n) in the state space has associated with it 
a certain probability between zero and one which is a function of time. 
In order to gain more insight into the stochastic process not a single 
realization of a battle has to be considered, but an average battle. 
Therefore the expected value of the force levels (i.e. averages) as a 
function of time and the variances in the force levels were investigated. 
The straight forward way to compute the expected force levels involves 


the knowledge of the probability distribution P(t,m,n). Then 


ue Ng 
E(M(t) > S__P(tama) (4.2) 
Mbp "bp 


n 
0 

N(t) > n P(t m,n) (as) 
"bp "bp 


are the expected values of the force levels as a function of time. 
There are some other ways to calculate the expected force levels, 

one of which will be stated here. Recall the “hybrid-analytical- 

numerical" method to get the state probabilities (equations 3.5 through 


th 


3.15). CLARK [4] has also shown that the i~ moment of, for example, 


the M-force level may be computed as 


= 


Ores 
E(M'(t) ) = 0f' } + = > 0) ; PE eaiebaien) (4.4) 


C4. 


with for 1SjSmp and 1$kSnq 


47 





j k 
em 2ochn 
Ske ay n=Q. vk eS), 
(7) 4 i O,n 
.* > iene (4.6) 


This again emphasizes the strong point of CLARK's solution. For a 
m,n 
j»k 
only once. Then with this information and relatively small computa- 


given set of battle parameters the coefficients C have to be computed 
tional effort not only the state probabilities but also the first and 
second moment of the force levels can be computed. This determines the 


variance in the force levels at the same time, e.g. 


Var(M(t) ) = E(Mo (t) ) - E(M(t))- E(M(t)) =. (4.7) 
On the other side, the weak point is that CLARK [4] considered only 
breakpoints equal Zero. 

Many authors have discussed one side's probability of winning or 
probability of winning conditioned on a certain number of survivors, 
which always eliminated the parameter time by integrating from time t=0 
to infinity. This may be legitimate to answer absolute (meaning time 
independent) questions about who will win, but for direct comparisons 
with the deterministic model, this author has the feeling that the best 
question to ask regarding a winner is: 

What is the probability of one side winning given the stochastic battle 
lasted as long as the deterministic one? 

The calculation of these probabilities gives another interesting 
probability, because given the time t = te (final time of the determin- 


istic battle) 
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P(battle has not yet finished) = 
1-P(M wins|t=t,)-P(N wins| t=t¢) (4.8) 


There has also been work done regarding the distribution of the 
time to finish a battle. But this is beyond the scope of this thesis 
(SPRINGALL [9]). 


SeeeOEFERENCES IN THE FORCE LEVELS 

The deterministic model, especially the force level equations (2.5) 
and (2.6), describe a process with continuous state parameters where, in 
reality, the possible states are integers. Quoting from LANCHESTER [7]: 

Since the forces actually consist of a finite number of finite units 
(instead of an infinite number of infinitesimal units) the end of the 
curve must show discontinuity, and break off abruptly when the last 
man 1s reached; the law based on averages evidently does not hold 
rigidly when the numbers become small. 

LANCHESTER suggested that his differential equations may be good 
approximations only as long as the force sizes are large. He also stated 
that the equations are based on averages, implying an underlying stochas- 
tic process. 

This shows that there must be a difference in the force levels which 


should become significant when the number of .combatants is small. This 


difference was called bias by CLARK [4] and TAYLOR [10]. It can be shown 


that 
GE(M(t)) = aE (N(t))+aB_ (t) Sa 
and 
dE(N(t)) _ 
MEL) = -bE(M(t)) +48, (t) wee ID) 
where 
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M9 Mg 
Be (t) = Mbp > P(t.msn)+ > nP(t.m,.n) (4.11) 
Mtl Nypt! 
and 
B(t) = M5 > P(t.m,.on)+ Y mP(t,m,ny,) . (2) 
Nopt! Mt! 


The bias terms B_(t) and BL (t) can be interpreted as the expected 
values of M(t) or N(t) conditioned on the fact that the battle has 


already ended at time t, for example, 
B(t)=E (M(t))[N(t)=n,, or M(t)=m, } (4.13) 


In other words, equation (4.9) says the expected casualty rate of the 
M-force is proportional to the expected number of survivors of the 
N-force given neither of the two forces has reached its breakpoint. 

Define the bias of the X-force as Ax(t) = E(M(t))-X(t) and the bias 
of the Y-force as Ay(t) = E(N(t))-Y¥(t). Then using equations (2.8) and 
(2.9) together with (4.9) and (4.10) it follows that 


d 
at 4* 


-adytaB (t) (4.14) 


d 
ap 4Y ~bAxtbB_(t) (4.15) 


with the initial conditions x(0)=) and Y(0)=). This has the solution 


t 
ax(t)=-Vab {18 (s) a/Bcosh(-Vab (t-s) 
0 
-B.(s)sinh(Vab (t-s))}ds, (4.16) 
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and 


le 
ay(t)= Yab_f 18,,(s) Yb/acosh( Vab(t-s) 
0 
-B.(s)sinh(-Vab(t-s))}ds . (4.17) 


Since for a fixed, nonegative argument z, the cosh(z) is always 
greater than the sinh(z), it is easy to visualize that in most of the 
cases both biases are positive, meaning the expected force levels of the 
Stochastic model are higher than the deterministic force levels. This 
has been shown by CLARK [4] and CRAIG [5] and confirmed by this author. 
In the rest of the cases the winner's bias is negative or close to zero 


and the loser's bias is positive. Two examples are given in Table 3. 
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FABER 3 


Examples for cases where one bias is positive and the other 


bias is negative. 


A. Y-force wins in a fight to the finish 


or (OC 40 de moO OG 

Xbp = 0 b = 0.04 

Yop = 9 

te = 15.581 Ax(ty) = 3.22 
dy (te) = -0.23 


B. Y-force wins in a fight with equal initial force levels, 


but different breakpoints. 


ne Le a = 0.08 

Xb p = 12 b = 0.08 

Ypyn = 8 

te = 2.85 Ax(ty) = 10.15 
Ay(te) = -0.10 


a2 





An interesting point has to be mentioned regarding case B of 
Table 3. Here the expected force level of the winner is 
smaller than the expected force level of the loser at the 
time a deterministic battle ends. 

To get a better feeling for the differences in the force 
levels, the next two figures show as an example a large spec- 
trum of force level behavior. The data and notation is 
described in Table 4. There are four battles with four dif- 
ferent breakpoints drawn as they evolve over time from zero 


to the time an equivalent deterministic battle ends. 
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X-force is always loser 


Y-force is always winner 
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corresponds to the battle with 


breakpoint force level: 
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Figure 15 - DETERMINISTIC AND EXFECTED FORCE LEVELS FOR 


THE X = FORCE 
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Figure 16 - DETERMINISTIC AND EXPECTED FORCE LEVELS FOR 


THE ¥Y - FORCE 
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CRAIG [5] has formulated hyootheses concerning the biases in the average 
force levels based on his work. His hypotheses were partially confirmed, 
but in some cases they have to be modified. Therefore, another similar 
set of hypotheses will be given and supported by Fig. 15 and Fig. 16, as 
well as some of the later figures. 
H 1) Given fixed initial force levels and attrition rate coefficients, 
as the breakpoint force levels increase, the numerical bias for the 
loser decreases. The biases for the winner do not show this monotone 
trend except for the case of symmetric parity. 
H 2) -Everything else constant, the bias of the loser increases with 
increasing initial force levels; this is also true in the symmetric 
parity case for both forces. 
H 3) Given the initial force level ratio is close to one at the time 
corresponding to the end of the deterministic battle, the bias of the 
loser is always larger than the bias of the winner. 


H 4) At the time corresponding to the end of a deterministic battle, 
the biases become larger as the forces come closer to parity. 


H5) The biases at times corresponding to less than one half the 
duration of the deterministic battle are negligible. 


The case of symmetric parity, i.e. equal initial force levels, 
breakpoints and equal attrition rate coefficients seems to be kind of a 
"limiting" case. For example, at parity the biases of both forces behave 
in the same manner and are equal. In Fig. 17 the biases at the time a 
deterministic battle ends as a function of the initial force levels and 
as a function of the breakpoints are presented. It is also another 


verification for the hypotheses H 1 and H 2. 
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Figure 17 - BIASES IN SYMMETRIC PARITY 
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Figure ig - BIASES WITH DIFFERENr INI1rIAL FORCE LEVELS 
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Figure 19 - BIASES WITH DIFFERENT aTTRITION RATE 
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Considering the changes in magnitude of the biases for battles like 
in Fig. 15, 16 and 17, CRAIG [5] came to the conclusion that when the 
forces are closer to parity, the biases at the deterministic battle's 
end increase. Several battles were fought starting with symmetric parity 
and then varying the force levels and the attrition rate coefficients. 
Sample results are shown in Fig. 18 and Fig. 19. In Fig. 18, the initial 
force levels were changed giving the Y-force ten and twenty percent 
higher initial force levels. The biases at the end of the deterministic 
battle are plotted as a function of the breakpoint force level ratio 
r= Xpp/ Xo: 

Fig. 18 supports CRAIG's hypothesis, as does Fig. 19. Here not the 
initial force levels but the attrition rate coefficients were changed in 
order to deviate from symmetric parity. The last way to deviate from 
symmetric parity is a case where the deterministic model gives equal 
answers for different battles. The data and the results are shown in 


table 5. 
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TABLE 5 


Equal initial force level battle with non equal breakpoints 


“7 alk aes ae a= b= f = 0.08 
Xb p ‘bp Ax(te) Ay(te) 
12 Zz 1.08 ey Ps: 
ip 9 Oey 1 D503 
UZ 6 O25 -0.10 
12 3 0.68 07258 

9 9 | Ao2 peor 
9 6 1.14 0.93 
9 3 1.06 0.75 
6 6 ili) Zeal 
6 3 1.67 ZZ 
3 3 2.84 2.84 
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The X-force level at time te 1S the same as the Y-force level. Thus, 
the expected force level for the winner (i.e. the Y-force) is smaller than 
the expected force level of the loser (the X-force), because ax(te) is 
always larger than Ay(te). This is easy to see when one remembers the 
way the expected force level is computed (equation (4.3)). Since the 
Y-force has a lower breakpoint, there are states (m,n) possible where 
Ny p<ricm, Apparently these states have a nonzero probability associated 
with them, which decreases the expected force level below the expected 
value for the X-force. This might be a starting point for further 


studies. 


D. VARIABILITY IN THE FORCE LEVELS 

Naturally in the deterministic case there does not exist any variabil- 
ity in the force levels. On the other hand, for the stochastic model, the 
variance in the force levels as a function of time is a measure of the 
dispersion of the number of survivors about their mean value. 

CLARK [4] has hypothesized two different types of behavior for the 
variance in the force levels, which are shown in Fig. 20 for the data 
presented in Table 2. The first type of behavior is that of the variance 
for the N-force, i.e. the variance increases monotonely as a function of 
time and is asymptotic to a limiting value. It was found that this type 
of behavior occurs when the side is going to win and for the case of 
symmetric parity. The second type of behavior shown 1s the variance of 
the M-force, Var(M(t)), as a function of time. This increases to a 
maximum value then decreases asymptotically to a limiting value. This 


type of behavior is associated with the loser of the battle. 
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The variance in the force levels is a function of the initial and 
breakpoint force levels, the battle time and the attrition rate 
coefficients. Unfortunately one does not know what this dependence is. 
Based on many numerical results a set of hypotheses will be stated and 
the next figures will support them. 

H 1) Given fixed initial force levels and attrition rate coefficients, 
as the breakpoint force levels increase the variance of the force levels 
decrease. 

H 2) Everything else constant the variance of the loser's force level 
increases with increasing initial force levels. This is also true for 
both variances in the case of symmetric parity. 

H 3) Given the initial force ratio is close to one at the time 
corresponding to the end of the deterministic battle the variance 

of the loser's force level is always smaller than the winner's variance. 
H 4) At the time corresponding to the end of the deterministic battle 
the variance in the loser's force level increases as the forces come 
closer to parity. This trend is not true for the variance of the win- 
ner, except for the case of symmetric parity. 

It was the intention of this author that the set of cases to illus- 
trate the hypotheses are the same as in the illustrations (Fig. 17, 18 
and 19) of the hypotheses about biases. So, Fig. 21 shows the force 
level variances for the different initial force levels and different 
breakpoints for the case of symmetric parity. In Fig. 22 the variances 
for a battle with equal attrition rate coefficients but varying initial 
force levels show that H 4 is only true for the loser. This point is 


emphasized by Fig. 23, where with constant and equal initial force levels 


the attrition rate coefficients were varied. 
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Figure 20 - rORCE LEVEL VARIANCES OVER TIME 
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Figure 21 - VARIANCES IN SYMMETRIC PARITY 
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Figure 22 = vARIANCES WITH DIFFERENT INITIAL FORCE LEVELS 
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Figure 23 - vARIANCES WITH DIFFERENT ATYRITION RATE 
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V. COMPUTATIONAL ASPECTS 


A. USING THE NUMERICAL SOLUTION 

For most of the numerical work, the state probabilities have been 
obtained using the fourth-order RUNGE-KUTTA method. The accuracy of the 
results was increased by substituting the available analytical results 
for the Regions I, II and III. 

The disadvantage is that this solution method needs a lot of CPU-time. 
For the battle in Example 1, with the data given in Table 1, the calcula- 
tion of the state probabilities, expected force levels and variances led 
to a CPU-time on an IBM-360 computer of almost 90 minutes with a time 
Step size of 0.05 minutes. That was the reason why, in the examples the 
attrition rate coefficients were increased by one magnitude, which brought 
the computer usage down to a CPU-time of around 12 minutes with the same 
time step size of 0.05 minutes. The sum of the state probabilities at 
every time step was used as a measure of accuracy. Surprisingly, its 
deviation from 1.0 was always in the fourth or fifth decimal, which proves 
the robustness of the RUNGE-KUTTA method. Also, the fact that for the 
Regions I, II and III analytical solutions were substituted at each time 
Step did not change the final outcome considerably. It was found that 
without the analytical partial solutions the sum of probabilities tended 
to be slightly higher, changing the sum of probabilities from, for 
example, 0.99998 to 1.00003. 
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B. THE EQUAL ATTRITION RATE COEFFICIENT SOLUTION 

After the development of the equal attrition rate coefficient solution 
(EARCS) outlined in Ch. III.C, it was implemented on an IBM-360 comput- 
er. The CPU-time for the calculation of the state probabilities, expected 
force levels and variances for a given point in time always stayed below 
20 seconds, which emphasizes its computational advantage. 

Further research showed two shortcomings of the EARCS, which are 
easily overlooked. The coefficients C(m,n) vary over a wide range start- 
ing at 1.0 and depending on the initial force levels. 


For example for a battle with 


My = 20 Mp = 15 
No = 40 Mbp = 20 
f= 0.08 te = 1.60995 


C(20,40) = 1.0 but (16,21) = 4.9297077276x10°°. 


Also, the binomial coefficients in equations (3.23) and 3.24) vary over a 
wide range starting at 1.0 to 1352078 for the above example. This shows 
that the capacity the computer, in terms of the number of significant 
digits, 1s able to carry limits the implementation of this solution in 
the present form. 

Another reason why this solution is not the end of the numerical 
problems is the fact that in the given form, equations (3.23) and (3.24) 
require computation of an alternate sum consisting of terms whose factors 
are relatively large binomial coefficients and small numbers between zero 
and one. This produces truncation errors, which yield nonsensical results 


like negative variances and sums of probabilities greater than one. 
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VI. CONCLUSIONS 


Throughout this thesis, only LANCHESTER-type square-law attrition 
with fixed initial force levels and fixed attrition rate coefficients 
has been considered for a deterministic and a stochastic version. The 
Stochastic version required much more computational effort. So, given 
the need for an analytical model as opposed to the use of simulation, 
there is not much to gain from the application of a stochastic model 
when the force levels are large and the forces are not near parity. In 
these cases the deterministic version essentially produces the same 
results, at least in qualitative terms. For smaller force levels or 
forces near parity, the stochastic version may be helpful to get more 
information about the dynamics of combat. 

Also, there exist one case where the deterministic version cannot 
differentiate between several different battles. That is the case of 
equal initial forces, equal attrition rate coefficients, but different 
breakpoint force levels. This also rectifies the further development of 
the equal-attrition-rate-coefficient-solution (EARCS). 

Using an analytical model, the discussed way of introducing random- 
ness into the model seems not to be very enlightening. Therefore, it 15s 
Suggested that another way to include random effects should be explored. 
The author's opinion is that working with attrition rate coefficients 
which are random variables seems more promissing to gain insight into 


the dynamics of combat. Further down the line there should be some 
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consideration on the usage of combinations of the possibilities to 


include random effects as given in Chapter II.B. 
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APPENDIX A 


THIS PROGRAM COMPUTES THE STATE PROBABILITIES USING THE FOURTH 
ORDER RUNGE-KUTTA METHOD. PLOTS OF THE EXPECTED FORCE LEVELS AND 


VARIANCES AS WELL AS 3-O-PLOTS FOR THE STATE PROBABILITIES ARE 
GBTAINED. 


DIMENSION Pi (33,33), Pe (33,33), 13D (20) 
DIMENSION XMT(SO00), XNT (SOO), QT(500), VM(S00), VN(S00), 


1 VMB(551), VNB(S51), VQ(551), VARX(S51), VARY (S51), TIM (S51) 


DIMENSION EST (313) ,0(313),F (2) ,SIZE (2) ,KX (100) , KY (100) ,WK (41,41, 3) 
i » PC (u1,41) ,CST (41) ,CETERM (41) 
LOGICAL™1 ION(Y1,41) 
DIMENSION DETERM (313) 

REAL™8 TTL (i2) Silex" 7A 
REAL Ki, Ke, K3, K4 

CALL ERRSET (208, 600,-1, 1) 

F (1) <0. 

F (2) =O. 

LINES=0 

SIZE (1) #6. 

SIZE (2) =8. 

QFLAG = 0. 

NKXY=100 

READ(S5,21) MBP, MO, NBP, NO 
READ(S,22) AA, BB 

READ (S,22) H, FTIM 

READ (5,21) N3D 

READ(S,22) EPSQ, EPSPTT 
T30(1) = 9999. 

IF (N3D .LT. 1) GO TO 25 
READ(S,22) (T3001), 1=#1,N3D) 


2S CONTINUE 


N30=NUMBER OF 3D-PLOTS, H=SIZE OF TIME STEP 
FTIME*FINALTIME, EPSQ,EPSPTT ARE ZERO LEVELS 
T3D(1) = FRACTION OF FTIME WHEN TO PLOT 3D 


16 
21 
2e 
28 
29 
30 


31 


DO 10 1#1,N30 

730 (1) sFTIM™T3D (1) 

CONTINUE 

FORMAT (1615) 

FORMAT (8F 10.5) 

WRITE (6,28) MBP. MO, NBP, NO 

FORMAT (// SX, ‘MS &N SS‘, SX, YI6 /) 
WRITE (6,29) AA, BB 

FORMAT (/ 5X, ‘A & Bt, 9X, 2F10.3 7) 

WRITE (6,30) 4H, FTIM 

FORMAT t/ SX, “INITIAL H & FINAL TIME (= LOOPS » H)*, SX, 2F10.3/) 
WRITE (6,31) N3D 

FORMAT (// SX, ‘s OF 3-D PLOTS*, 2x, 16 /) 
IF (N3D .LT. 1) GO TO 36 
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c¥y 


32 
36 


47 


48 


50 


52 


Ss 


WRITE (6,24) EPSO,EPSPTT 
PORMAT(, “5x, ERSO=°,F10.5.5X. “EPSPIT=*,F10.5) 
WRITE (6,32) ({T30(1), J=1,N3D) 
FORMAT (/ SX, ‘AT TIME’, 2X, 10F10.3 7) 
CONT INUE 

MO1 = MO + 1 

NOi = NO + 1] 

NS = NOI - NBP 

MS = M01 - MBP 

MSL = MS + } 

NSL = NS) +1 

BMG = BB = MOQ 

ANG = AR = NO 

QANB = ANO / BB 

QBMA = BMG / AA 

8AT = -{(BMO + ANO) 

00 4? M1 = 1,M5 

YM(M1) = MBP + M1 - 1 
CONTINUE 

DO 4¥8B Ni = 1,NS 

VNI(N1) = NBP + NI - 1 
CONTINUE 

L = 0 

N3 = 1 

TPLOT = 130 {N3) 

TIME = OQ. 

CONTINUE 

IF (TIME .GT. FTIM) GO TO 210 
CONT I NUE 

TIME = Hw L 

et 

L8 20 

TIM{L) = TIME 

IF ( & .EQ. 2) GO TO 53 

IF (TIME .LT. TPLOT) GO TO 60 


SET NEXT PLOT TIME 


N3 = N3 + J 

TPLOT = 130 (N3) 

IF (N3.GT.N3D) TPLOT=9999. 
CONTINUE 


L8 = &£ WHEN TIME TO PLOT 


L8 =#L 

E8T = EXP(BB x TIME) - 1. 
EAT = EXP{AA x TIME) - 1. 
AE = QANB = EBT 

BE = QBMA = EAT 
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101 
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104 
105 
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108 
109 
110 
111 
lie 
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116 
11? 
118 
119 
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121 
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123 
124 
125 
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128 
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132 
133 
134 
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136 
137 
138 
139 
140 
141 
142 
143 
gy 
145 
146 
147 
148 
149 
150 


a | 


i ae J 


Oo 


a 


60 


502 


BAE = EXP(BAT x 
FJ = 1. 
JJ = 1 


NO VALUE OF J = 


FK = 1. 
KK = | 
QF Ke}. 
QFJ=1. 


NO VALUE OF K 


CONTINUE 


IF 


{(L .EQ. 1) G 


QQ=0. 

00 4OO N1 = 1,N 
N = NOI - Ni 

NL = NSL - Ni 
NL1 = NL + 1 
ANJ = AR =» N 


TIME) 


ZERO AT MO,NO 


ZERO AT MO,NO 


0 TO 500 


S 


00 300 M1 = 1,MS 
M = M01 - M1 

ML = MSL - M1 
ML1 = ML + 1} 
BMI = BB = M 
ABA = ANJ + BMI 


IF 
IF 
IF 
IF 


K1 2 ANJ*=Pe (ML1,NL) 
PT = ANJI=O. Sx 


(M .EQ. MQ) 
(4M .EQ. MBP) 
(N .EQ. NO) 
(N .EQ. NBP) 


GO 10 502 
GO TO SOY 
GO 70 506 
GO 10 507 


DEFAULTS TO ALL INSIDE POINTS 


1 Pe (ML,NL1)) 
K2 = PT - ABA x 


K3 = PT - ABA = (P2(ML,NL) 
KY = ANJ»«P1 (ML1,NL) 


P1I(ML,NL) = Pel 
G6 TO 200 


IF 
IF 


(N .EQ. NO) 
(N .EQ. NBP) 


M=MO, N 


+ BMI»P2 (ML,NL1) ~- ABAxPe (ML.NL) 
(P1(ML1,NL)+P2(ML1,NL)) + BM]=*O.S =» (P1(ML,NLI) 


(P2(ML,NL) * HxO.5™K1) 


ML,NL) + (H/6.0) x 


GO TO 503 
GO 10 507 


—=N0, NBP 


K1 = BMI™P2(ML,NL1) - ABAxPe (ML, NL) 


75 


+ HxO.5=K2) 
+ BMI=P1(ML,NL1) - ABAx(P2(ML,NL) + HxK3) 


(K1it2.0«Ke + 2.0»K3 + KY) 


+ 





151 
152 
153 
154 
155 
156 
157 
158 
159 
160 
161 
162 
163 
164 
185 
166 
167 
168 
169 
170 
171 
172 
173 
174 
175 
176 
177 
178 
179 
180 
181 
182 
183 
184 
185 
186 
187 
188 
189 
190 
191 
192 
193 
194 
195 
196 
197 
198 
199 
200 


O00 On-0O0 aan 


on 


oon 


503 


504 


505 


506 


507 


PT = BMIéO.5S = (P1(ML,NL1) + P2(ML,NLI)) 

Ke = PT - ABA = (P2(ML,NL) * HxK1~0.5) 

K3 = PT - ABA = (P2(ML,NL) + H¥Ke=x0.5) 

KY = BMI = PI(ML,NL1) - ABA = (P2(ML,NL) + HxK3) 

P1(ML,NL) = P2(ML,NL) + (H/6.) »« (K1it2.0xK2 + 2.0«K3 + KY) 
}FetLS .EG. L) GA 16 60! 

GO TAG 200 


M=M40, N=NO 


CONTINUE 

Pi(ML,NL) = EXP(-ABA = TIME) 

WRITE (6,102) L, N, ML,NL, Pi (ML,NL) 
GO 10 200 

IF (N .£Q. NBP) GO TO 505 


M=MBP, N-= NBP 


Ki = ANJ » P2(ML1,NL) 

K2 = ANI = 0.5 = (P1(ML1,NL) + P2(ML1I,NL)) 

K3 = Ke 

KY = ANS = Pi (ML1,NL) 

Pi (ML,NL) = P2(ML,NL) + (H/6.0) = (Ki + 2.0xK2 + 2.0xK3 + KY) 
GQ TG 200 


M=MBP, N=NBP 


CONTINUE 
P1(ML,NL) = 0.0 
GO 170 200 
N=NO, M-=MBP, MO 
Ki = ANJ«P2(ML1,NL) =- ABAxP2 (ML,NL) 


PT = ANS =» 0.5 = (PI(ML1,NL) + P2(MLI,NL)) 

K2 = PT - ABA = (P2(ML,NL) + HxO.SK1) 

K3 = PT - ABA = (P2(ML,NL) + HxO.5K2) 

KY = ANS » P1(ML1,NL) - ABA = (P2(ML,NL) + HxK3) 

PI(ML,NL) = PO(ML,NL) + (H/6.0) = (K1 + 2.0xK2 + 2.0™K3 + KY) 
IF (LB .EQ. L) GO TO 602 

GO TO 200 


N=NBP, & N=NBP, M=MO 


CONTINUE 

Ki = BM] =» PO (ML,NL1) 

Ke = BM] =» 0.5 = (PI(ML,NL1) + P2(ML,NL1)) 
K3 = Ke 

KY = BM] = PI(ML,NL1) 


76 





201 
202 
203 
204 
205 
206 
207 
208 
209 
210 
elt 
ele 
213 
214 
215 
216 
217 
218 
219 
220 
ee! 
eee 
223 
224 
225 
226 
2e7 
226 
229 
230 
231 
232 
233 
234 
235 
236 
237 
238 
239 
240 
2ul 
ede 
2u3 
euy 
2us 
246 
2h? 
248 
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601 


250 


259 
260 


602 


265 
270 
275 


700 


701 


PIUML. NG "Fem. NL) + (H/6.0) = (KI * 2.0“K2 + 2.0%K3 + KY) 
GO 10 200 


ALSO M=M0, N-=NBP, NO 


CONTINUE 
IF C(BE.LT.1.) .AND. (KK.GT.9)) GO TO 255 


COMPUTE PTT IF BE > O & KK < 10 


PTT=BAE 

06 250 IND=1,KK 
XIND=INO 
PTT=PTT*BE/XIND 
CONT I NUE 

GO TO 260 
PTT=0. 

CONT INUE 

KK = KK + 1 

IF (CABS (P1(ML,NL) -PTT)).LE.EPSPTT) GO TO 199 
GO TO 700 


M=MBP, N=NO 


CONTINUE 
IFC (AE.LT.1.) .AND. (JJ.GT.9)) GO TO 270 


COMPUTE PTT IF AE > O & JJ < 10 


PTT=BAE 

00 265 IND=1, JJ 
XIND=IND 
PTT=PTT»AE/XIND 
CONTINUE 

GO TO 275 
PTT=O. 

CONTINUE 

JJ 2 JS *+ 1 

IF (CABS (P1 €ML,NL) -PTT)).LE.EPSPTT) GO TO 199 
CONTINUE 


TIME REOUCED 1/2 


Lw2zet - 1 

WRITE (6,701) TIME, M, N. H 

FORMAT (// SX, ‘H VALUE 15 REDUCED BY HALF AT : TIME = °*, F8.3, 
ene ° 15,  Ne-= *, 15, ° FROM *. F8.3 //) 

H= 0.5 * H 

WRITE (6,102) LL.M.N.KK, PTT, PI CML.NLI,QOFK, QFJ,00,QFLAG 
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gol 
252 
253 
254 
255 
256 
257 
258 
259 
260 
261 
262 
263 
264 
265 
266 
267 
268 
269 
270 
27] 
27e 
273 
e7y 
275 
276 
277 
278 
279 
280 
281 
282 
283 
28y 
285 
286 
287 
288 
289 
290 
Zo! 
292 
293 
egy 
e95 
296 
297 
298 
299 
300 


mo 


om 


Om 


102 


199 
200 


300 
400 


401 


¥10 
W15 


FopmMet (7 2X, 41S, 7Fi2.5) 
QFLAG=0. 

G6 TO 52 

PTT=P1 (ML,NL) 

CONTINUE 

QQ=Q0+P1 (ML,NL) 

CONTINUE 

CONTINUE 

Q01=ABS (QQ-1.) 

IF (QQ1.LE.EPSQ) GO TG 401 
QFLAG=1. 

G6 T6 700 

CONT INUE 


COMPUTE VARX(T), VARY (1) ,N3SAR (1) , MBAR (T} ,O (7) 


SmM8 = QO. 

SN8 = QO. 

$0 = 0. 

SM2 = 0. 

SNe = 0. 

06 41S MI = 1,M5 

RM = VM (M1) 

RM2 = RMxRM 

06 Y¥10 Ni = 1,NS 

RN = VNINI) 

RN2 = RNxRN 

PT = P1(M1,N1) 
PO(M1,N1) = PT 

SM8 = SM8 + RMxPT 
SN8 = SN8 + RNxPT 

SQ = SQ + PT 

SM2 = SM2 + RM2 = PT 
SN2 = SN2 + RN2 = PT 
CONTINUE 

CONTINUE 

VM8(L) = SMB 

VNB (L) = SNB 

VO(L) = SQ 

VARX(L) = SM2 - SMBxSMB 
VARY (L) SNe - SNBxSNB 


RETURN TO MAIN LOGP (50) IF NOT TIME TO PRINT 
IF (N30 .LT. 1) GO TO 430 

IF (L8.NE.L) GO TO 430 

IF (L8.EQ.2) GO TO 430 


ROJUST Xst, VECTORS FOR PLOT 
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301 D0 15 K=1,M5 


302 1S EST (K) =K-1 

803 DO 16 K=1,NS 

304 DE TERM (K} =K-1 

305 16 CONTINUE 

306 DO 115 I#1,41 

307 CST) #] 

308 CETERM (1) =] 

309 DO 115 J=1,41 

310 PC {I,J} #0.0 

311 115 CONTINUE 

312 MP=41 

313 NP=41 

314 ALP=#15. 

315 BETA=30. 

316 PMAX=P1 (1,1) 

317 00 116 1=1,NS 

318 WRITE (6,199) 1,P1(1 ,1)},P101,1 ) 
319 199 FORMAT(* *,15,10X,2F15.5) 
320 116 CONTINUE 

321 DO 18 121,MS 

322 DO 18 J=1,NS 

323 IF (P1(1,J).GT.PMAX) PMAX=P1 (I,J) 
324 18 CONTINUE 

325 CONST=8.0/PMAX 

326 DO 19 I=1,MS 

927 N=] +MBP 

328 DO 19 J=1,NS 

329 M=J+NBP 

330 P1 (I,J) =CONST™P1 (1, J) 
931 PC(N,M)=P1 (1, J) 

332 19 CONTINUE 

333 CALL PLT301 (CST,MP,CETERM,NP,PC,ALP,BETA,F,TIL. SIZE. WK, ION, KX,KY,N 
334 IKXY,LINES) 

335 6 

336 430 CONTINUE 

337 & 

338 c 

339 c PLOT HERE 

340 c RETURN TO MAIN LOGP 
9u1 fe 

942 Go To 50 

343 S00 CONTINUE 

3uy DO SO1 1 = 1,MS 

3u5 DO SO1 J = 1,NS 

346 P1{1,J) = 0.0 

3u7 S01 CONTINUE 

3u8 P1(MS,NS) = 1.0 

3u9 GO TO 401 

350 c 
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351 
352 
353 
354 
355 
356 
357 
358 
359 
360 
361 
362 
363 
364 
365 
366 
36? 
368 
369 
370 
371 
372 


FINAL TIME REACHED 


210 CONTINUE 
WRITE (6,211) TIME, FTIM, L 
211 FORMAT(// 5X, “COMPUTED TIME ‘*, F10.3, 5X, ‘INPUT FINAL TIME °, 
1 F10.3, 5X, ‘** OF LOOPS TO REACH FINAL TIME *, I5 //) 
00 215 I =1,L 
WRITE (6,213) VMB(I), VNB(I), VO(1), VARX (I), VARY (I), TIM({T) 
213 FORMAT (2X, 12F10.5) 
215 CONTINUE 
650 FORMAT(*! ) 
WRITE (6,650) 
CALL PLOTP (TIM, VMB,L, 0) 
WRITE (6,650) 
CALL PLOTP{TIM, VARX,L,0) 
WRITE (6,650) 
CALL PLOTP (TIM, VNB,L, 0) 
WRITE (6,650) 
CALL PLOTPI{TIM, VARY,L,0) 
STOP 
DEBUG SUBCHK 
END 
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oor ouvwlwnhrv = 


onononn 


APPENDIX B 


THIS PROGRAM CALCULATES THE STATE PROBABILITIES FOR THE EQUAL 
ATTRITION RATE COEFFICIENT SOLUTION {ERRCS). 


A 


3-D PLOT 1S PRODUCED USING THE VERSATEC PLOTTER. 


IMPLICIT REALx8 {A-H,0-Z) 
CALL ERRSET (208,256,10,1) 
DIMENSION C (50,50) ,4(50) ,N(50) ,F (50) ,PTMN (SO, 50) 
REAL™Y SIZE (2),FL (2), RK (41,41,3),X (41), 7041), P41, 4d) 
DIMENSION KX (100) ,KY (100) 
LOGICAL™! ION (41,41) 
REAL*8 TTL (12) /12x* 4 
READ (5,100) M0, MBP,NO.NBP 
READ (S,101)A 
READ (S, 102) TIME 
READ (5,103) (F (1) ,1=1,41) 


F(I)=1/1-FACTORIAL, DONE TO SPEED UP THE PROGRAM 


802 


10 


20 


21 


WRITE (6,802) 40,MBP,NO,N8P,A,TIME 
FORMAT(* ‘*,415,2F10.5) 
MD=MO-MBP 

ND=NO-NBP 

MO1=MO+1 

NO1=NO+1 

MBP1=MBP+1 

NBP1=NBP+! 

MBP2=MBP +2 

NBP2=NBPte2 

RMO=MG 

RNO=NO 

08 10 1=1,41 

DO 10 J=1,41 
C(1,J)=0.0 
PTMN (1, J) =0.0 
P(I,J)=0.0 
CONTINUE 

C (M01,N01) =1.0 

D8 20 1=2,mD 

J=MBP +] 

C(J,NO1)= RNOx™ (MD-1+1) 
CONTINUE 

DO 21 I=2,ND 

J=NBP+] 

C(MO1,J)= RMOx™~ (ND-141) 
CONTINUE 

DO 22 1=2,MD0 

MC=M01-1 

MM=MO1-1+1 

MPLUS=MM+1 
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ee 


C END 


25 
100 
101 
102 
103 


C FOR 


302 


301 


30 


DO 22 J=2,ND 

NC=NO1I-J 

NN=NO1-J*! 

NPLUS=NNF1 

C(MM,NN) =NC= C(MPLUS,NNJ +MC = C(MM,NPLUS) 
CONTINUE 


OF COEFFICIENT CALCULATION 


Fi2-AxTIME 

DO 25 1=MBP2,M0! 

DO 25 J=NBP2,NO0! 
]J=l+J 

K=MO+*NO-IJ +2 

Fe=(1.0-DEXP (F1)) «kK 
KK=K+1! 

IJ=]J-2 

F3=DEXP (Fixl J) 

PTIMN (I,J) sF2xF3=C (1, J) «F (KK) 
CONTINUE 
FORMAT (415) 
FORMAT (F 10.5) 
FORMAT (F 10.5) 
FORMAT (017.11) 


NBP<N<NO 


DO 30 Il=!,ND 
NFORCE=NBP+tI 1 
J=MO0+NO-MBP-1-NFORCE 
FACT=MBP+!1.0+NFORCE 
SUM1= (1.0-DEXP (FlxFACTII/ ( AxFACT) 
SUMINT=SUM! 
DO 301 K=1,J 
FACT=FACT+1.0 
RODFAC= (-1.0) xxK 
FRACTN= (1.0-DEXP (FixFACTI) 7 ¢€ A»FACT) 
COMBT=1.0 
DO 302 KJ=1,K 
RKJ=KJ 
COMBI = {(J-RKJ+1.0) /RKJ 
COMBT=COMBT=COMBI 
CONTINUE 
SUMINT=SUMINT*ADDFAC=FRACTN=COMBT 
CONTINUE 
NMBP=J+! 
NBOUND=NFORCE*1! 
PFAC=RxNFORCE*F (NMBP) *C (MBP2, NBOUND) 
PTMN (MBP1,NBOUND) =PFRAC=SUMINT 
CONTINUE 
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101 
102 
103 
104 
105 
106 
107 
108 
109 
110 
111 
112 
113 
114 
115 
116 
117 
118 
119 
120 
121 
122 
123 
124 
125 
126 
127 
128 
129 
130 
131 
132 
133 
134 
135 
136 
137 
138 
139 
140 
141 
142 
143 
14y 
145 
1u6 
iu? 
148 
149 
150 


Cc 


C FOR MBP<M<MO0 


C 


c 


402 


401 


40 


801 
45 


BDO YO Il=1,MD 
MFORCE=MBP+t1 1 
J=40*+NO-NBP-1-MFORCE 
FACT=NBP+1.0+*MFORCE 
SUM1= (1.0-DEXP (FilxFACT)) /€ AxFACT) 
SUMINT=SUM1 
00 401 K=1,3 
FACT=FACT+1.0 
ADDFACs= (-1.0) «xk 
FRACTN= (1.0-DEXP (FixFACTI)/( AxFACT) 
COMBT=1.0 
DO YOe Kj=1,K 
RK J=KJ 
COMBI= (J-RKJ+1.0) /RKJ 
COMBT=COMBT=COMB! 
CONTINUE 
SUMINT=SUMINT+AODFACxFRACTN»COMBT 
CONTINUE 
MNBP=J+1 
MBOUND=MFORCE*! 
PFAC=AxMFORCE™F (MNBP) »C (MBOUND, NBPe) 
PTMN (MBOUND, NBP1) =PFAC™SUMINT 
CONT INUE 
DO YS 1=MBP1,401 
D6 YS J=NBP1,NO1 
K=]-1 
L=J-1 
WRITE (6,801) K,L.PTMN(I, J) 
Bontalt: P(l,”, 15,°,°.15,°) = *.D17.11) 
CONT INUE 


C DATA ADJUSTMENT FOR PLOT 


c 


= 


NROW=41 
NCOL=41 
NKXY=100 
LINES=0 
ALPHA=15. 
BETA=30. 

FL (1) =0.0 
FL (2) =0.0 
SIZE (1) =6.0 
SIZE (2) =8.0 


C SCALING 


Cc 


PMAX=PTMN (1,2) 
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151 
152 
153 
154 
155 
156 
157 
158 
159 
160 
161 
162 
163 
164 
165 
166 
167 
168 
169 
170 
171 
172 
173 
174 
175 
176 
177 
178 
179 
180 
181 
182 
183 
164 
185 
186 
187 
188 
189 
190 
191 
192 
193 
194 
195 
196 
197 
198 
199 
200 


Ooan 


Cc 


DO 50 I=! 
X (1) =] 
DO 50 J=! 
Y (J) =J 


»4)1 


41 


IF (PTMN (1, 3) .GT.PMAX) PMAX=PTHMN (I, J) 


50 CONTINUE 


CONST=8.0/PMAX 


DO 51 I=! 
DO 51 J=! 


91 CONTINUE 


41 
41 
P(I, J} = SNGL (CONST»PTMN (1, J)) 


IF OTHER FORCE LEVELS CHANGE PLOT ARGUMENTS 


CALL PLT301 (X, NROW, Y,NCOL,P,ALPHA, BETA,FL,TIL,SIZE,WK, 
ION, KX,KY,NKXY,LINES) 


EXPECTED VALUES AND VARIANCES 


EmM=0.0 
EMM=0.0 
VARM=0.0 


SPROB=0.0 
DO 60 I1=MBP1,401 


J=]I-1 
SPT=0.0 


DO 601 JJ=NBP1,NQ01 


SPT=SPT+ 
601 CONTINUE 


PTMN (11, JJ) 


EM=EM+] «SPT 
EMM=EMM+ I =I xSPT 
SPROB=SPROB+SPT 


60 CONTINUE 


VARM=EMM-EM™EM 


EN=0.0 
ENN=0.0 
VARN=0.0 


DO 70 JI=NBP1,NO01 


I=]J-1 
SPT=0.0 
DO 701 JJ 
SPT=SPT+ 
701 CONTINUE 


=MBP1,M01 


Finn, 11} 


EN=EN+1 SPT 
ENN=ENN+ II] SPT 


70 CONTINUE 


VARN=ENN-EN*EN 
WRITE (6,805) EM,EN, VARM, VARN 


805 FORMAT (* 


Pe EXE COTEDRVALUES: HM 
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N*, 


SOX, * VARIANCE M,N 


et A 





201 
20e 
203 
204 
205 


Meme Oem eOr. O17.11,151,017.11, 10X. 
WRITE (6,1100) SPROB 

1100 FORMAT (* SUM GF PROBABILITIES *,017.11) 
STOP 
END 
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